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We perform simulations of relativistic binary stars in post-Newtonian gravity to investigate their 
dynamical stability prior to merger against gravitational collapse in a tidal field. In general, our 
equations are only strictly accurate to first post-Newtonian order, but they recover full general rel- 
ativity for spherical, static stars. We study both corotational and irrotational binary configurations 
of identical stars in circular orbits. We adopt a soft, adiabatic equation of state with F = 1.4, for 
which the onset of instability occurs at a sufficiently small value of the compaction M/R that a 
post-Newtonian approximation is quite accurate. For such a soft equation of state there is no inner- 
most stable circular orbit, so that we can study arbitrarily close binaries. This choice still allows us 
to study all the qualitative features exhibited by any adiabatic equation of state regarding stability 
against gravitational collapse. We demonstrate that, independent of the internal stellar velocity 
profile, the tidal field from a binary companion stabilizes a star against gravitational collapse. 
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I. INTRODUCTION 

Binary neutron stars are known to exist and for some 
of these systems in our own galaxy (including PSR 
B1913-I-16 and B1534-I-12), general relativistic effects 
in the binary orbit have been measured to high pre- 
cision Interest in binary neutron stars has been 
stimulated by the prospect of future observations of ex- 
tragalactic systems by gravitational wave interferometers 
like LIGO |], VIRGO g, TAMA g and GEO §. Bi- 
nary neutron stars are among the most promising sources 
of gravitational waves for these detectors, and therefore 
it is important to predict theoretically the gravitational 
waveform emitted during the inspiral and the final coa- 
lescence of the two stars. Interest in these systems also 
arises on a more fundamental level, since the two-body 
problem is one of the outstanding unsolved problems in 
classical general relativity. 

Considerable effort has gone into understanding binary 
neutron stars. Most of this work has been performed 
within the framework of Newtonian and post-Newtonian 
gravity (see, e.g., for a review and list of references). 
General relativistic treatments are currently only in their 
infancy. Recently, Wilson, Mathews and Marronetti ||] 
(hereafter WMM) reported results obtained with a rela- 
tivistic hydrodynamics code. Their code assumed several 
simplifying physical and mathematical approximations. 
Their results suggest that the central densities of the 
stars increase as the stars approach each other and that 
massive neutron stars, stable in isolation, individually 
collapse to black holes prior to merger. WMM therefore 



find that in general relativity, the presence of a compan- 
ion star and its tidal field tend to destabilize the stars in a 
binary system. This conclusion is contrary to what is ex- 
pected from Newtonian , post-Newtonian [l0| - [l^, per 



turbative [ ]13[ and matched asymptotic expansion []14|, 15 
treatments of the problem. Constructing self-consistent, 
fully relativistic initial data for two neutron stars in a 
circular, quasi-equilibrium orbit does not show any ev- 
idence of this "crushing effect" either jl^]. Moreover, 
applying energy turning-point methods to sequences of 
these initial data suggests that inspiraling neutron star 
binaries are secularly stable all the way down to the in- 
nermost stable circular orbit |l7|| . To summarize, most 
researchers currently believe that the maximum allowed 
rest mass of neutron stars in close binaries is larger than 
in isolation, and that their central density is smaller than 
in isolation. If there exists any destabilizing, relativistic 
effect at high post-Newtonian order, then this effect is 
much smaller than the dominating stabilizing effect of 
the tidal field. 

However, to date, the only fully dynamical treatment 
of the problem in general relativity - that of WMM - 
reports a star-crushing effect. In this paper, we perform 
a new, fully dynamical simulation for binary stars in post- 
Newtonian gravity. We use a formalism in which (1) all 
first post-Newtonian terms are taken into account, and 
(2) sufficient nonlinearity is retained, so that spherical, 
static stars satisfy the fully general relativistic equations 
exactly. As explained in section II below, this formalism 
is very suitable for studying binary neutron stars. We 
study relativistic effects in binary stars with M/R ^ 1, 
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where M and R are typical values of the stellar mass and 
radius, so that a post-Newtonian treatment is completely 
adequate. 

By performing a fully dynamical calculation, we can re- 
lax various constraints assumed in previous treatments. 
For example, Wiseman assumed the stars to remain 
spherically symmetric | [T^ ] , Baumgarte et al. assumed 
the binary stars to be corotating, and Thorne as- 
sumed the stars' orbital separation to be much larger 
than the stars' radius. Here, we relax all these assump- 
tions and study tidally deformed stars, both corotational 
and irrotational, at arbitrarily small separations. We still 
find that the presence of the tidal field of a companion 
star tends to stabilize neutron stars against catastrophic 
collapse. 

To establish the stability of binary stars against col- 
lapse, we construct quasi-equilibrium initial data for 
identical binary neutron stars in a close, circular orbit. 
The idea is to show whether stars in a binary formed 
from the inspiral of objects which are stable in isolation 
remain stable at close separation. Our models have rest 
masses near the maximum allowed rest mass for spherical 
stars in isolation and thus provide the best candidates for 
collapse if the tidal field is destabilizing (stars with rest 
masses well below the maximum allowed value are un- 
ambiguously stable) . In order to demonstrate that these 
stars are dynamically stable, we need to locate the on- 
set of instability in the binary, and compare it with the 
onset of instability for isolated stars. Since the shift is 
fairly small, a very careful treatment with high numerical 
accuracy is necessary. We detail our method of locating 
the onset of instability in section IV. 

The paper is organized as follows: In section II, we 
present the post-Newtonian formalism adopted in this 
paper. We cahbrate our code in section III by locating 
the analytically known onset of radial instability of rela- 
tivistic spherical stars against gravitational collapse . 
In section IV, we study the dynamical stability against 
gravitational collapse of close binary stars, and briefly 
summarize our results in section V. 



II. FORMULATION 

In the usual post-Newtonian treatment, the fluid and 
field equations are derived by systematically expanding 
the Einstein equation and the relativistic hydrodynamic 
equations in powers of [ p^ . In this paper, we intro- 
duce a different approach. Since in a first order post- 
Newtonian approximation, the spatial metric may 
always be chosen conformally flat, we can derive post- 
Newtonian equations by starting with the relativistic 
equations written previously in the conformally flat ap- 
proximation We then neglect some of the second 
and higher order post-Newtonian terms, but retain suf- 
ficient nonlinearity so that this formalism recovers full 
general relativity for some limiting regimes of interest in 



this paper. 

We write the spatial metric in the form 

li] = ^^li] = 'ip'^S^j, (2.1) 

so that the line element becomes 

ds^ — Qfii^dx^dx'^ 

= {~a^ + f3kP^)df + 2(3dx'dt + i^^S.jdx^dx' . (2.2) 

Here, a, [3i and ip are the lapse function, the shift vector 
and the conformal factor and we adopt geometrized units 
in which c = 1 = G. We also adopt cartesian coordinates 
x^ — {x,y, z), so that the covariant derivative Vfe associ- 
ated with jij = 6ij conveniently reduces to the ordinary 
partial derivative d/dx^ . 

We employ a perfect fluid stress-energy tensor 

T^" = p (^1 + £ + u^u'' + (2.3) 

where p, e, P, and denote the rest mass density, spe- 
cific internal energy, pressure, and the fluid four velocity, 
respectively. For initial data, we assume a constant en- 
tropy configuration with a polytropic pressure law 



P = Kp^ 



(2.4) 



where K and F = 1 + 1/ri are constant and n is the 
polytropic index. For the evolution of the matter we 
assume the adiabatic relation 



P = (F - l)pe. 
The continuity equation is 
9p* d{p^v^) 



dt 



= 



(2.5) 



(2.6) 



where p^^ = pau^ip^ and — /u^ . 
The relativistic Euler equation is 

d{p^u,) d{p^UiV^) 6d ~o 



dt 



dx^ 



~ r,] , '2p*UkUk , 

■^P*u,P\^+ (2.7) 

and the energy equation is 

9e* d{e^,v^) 



dt dxi 



= 0, 



where 



Uj = {1 + Te)uj, 
u" = (1 + Fe)u", 



e* — (p*e*)^^'"; e* = e(au'V°) 



o„/.6Nr-i 



(2.8) 



(2.9) 
(2.10) 
(2.11) 
(2.12) 
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Note that /?' = ip~'^f3i in the conformal flat approxima- 
tion. From the normahzation condition u''Ui = —1, we 
find 



[au 



0\2_ 



1 



1 



UjUj 



1 + r 



o./,6^r-i 



(2.13) 



In our nmncrical simulation, we use p*, Ui, and as the 
independent variables that are determined by the hydro- 
dynamical equations. 

Equations for ip, f3^ and a can be found from the 
Hamiltonian constraint, the momentum constraint, and 
the maximal slicing condition tvK = 0, where tvK is the 
trace of the extrinsic curvature Kij (details can be found, 
for example, in 

Assuming maximal slicing for all times, we have 
dttvK — 0, and can use the trace of the evolution equa- 
tion for Kij to find 

(2.14) 

Here E = p{l + Te){au°)'^ - P, Sij = Tij, and A is the 
flat space Laplacian. 

Since we assume the spatial metric to remain confor- 
mally flat for all times, the trace free part of the time 
evolution equation for jij has to vanish, which yields 



(2.15) 



This equation shows that the extrinsic curvature no 
longer represents independent dynamical degrees of free- 
dom (i.e., it may no longer exactly satisfy its fully rela- 
tivistic evolution equation). Inserting this into the mo- 
mentum constraint, {tp^K^j)^i = S^Jj-ip^, where J; = 
p^,Uiip~^ , we find 



A/3* 
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In 



(^)l i^^ 



Finally, the Hamiltonian constraint yields 



(2.16) 



(2.17) 



For the post-Newtonian point of view assumed here, 
a conformally flat spatial metric takes into account all 
Newtonian and first post-Newtonian terms, and differ- 
ences from a general, fully nonlinear metric appear at 
second post-Newtonian order pO| , ^ . We can therefore 
use the above equations, which assume a conformally 
flat metric, as a starting point for a post-Newtonian 
approximation. Wc simplify the problem by neglect- 
ing other terms of second or higher post-Newtonian or- 
der. In particular, we will neglect the nonlinear terms 



^(3\.^-2S,j(3\i/3), and 
KijK'^^tfj^ /8. Note that for static, spherically symmetric 
spacetimes, these terms vanish identically, so that we still 
recover full general relativity for these spacetimes. 

Adopting this approximation, the field equations re- 
duce to 



A(ai/;) = 2TTai^^(^E + 2S^jS'^tp-^^ 
A/3* 



/3"' = 167ra J,; 



Ai; = -2iTil>'^E EE 47r5^. 



We decompose the equation for using 



/3* = 45, 



so that Bi and x satisfy 



ABi = AivaJi , 
Ax — —AiraJiX^ 



(2.18) 

(2.19) 
(2.20) 

(2.21) 



(2.22) 
(2.23) 



To summarize, we have reduced Einstein's equations to 
six elliptic equations for the six functions aip, tp, Bi and 
X- We solve these equations together with the boundary 
conditions 



V-- 1- - / S^dV + 0{r-'), 



(2.24) 
(2.25) 



B^ 



aJxxdV — ^ / aJxljdV 



+0(r-4), (2.26) 

X f y f 

= J o^JyxdV - ^ y aJyydV 

+0(r-^), (2.27) 
B^^ j aJ^zdV + 0{r-^), (2.28) 

j ajydV + 0{r-'^), (2.29) 

where dV is the coordinate volume element. Note that 
having removed s ome of th e non linear terms from the el- 
liptic equations ( 2.18] ) to ( 2.23 ), their right hand sides 
now have compact support. This further simplifies the 
computations, since in imposing the boundary conditions 
at a finit e sep arat ion, we do not truncate any source 
terms in ( 2.18| ) to (2.23) that extend to infinity. 

Having neglected some second post-Newtonian terms, 
our formalism is strictly only first-order post-Newtonian. 
Note, however, that we have only truncated some of the 
non-linear terms in the field equations, pieces of which, 
loosely speaking, can be associated with dynamical fea- 
tures of the gravitational fields. In particular, we still 
solve the fully relativistic hydrodynamic equations. We 
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therefore retain many of the nonlinear features of full gen- 
eral relativity, and expect that this formalism provides 
an excellent approximation in several limiting regimes of 
interest here. 

For example, for static, spherically symmetric stars, 
we recover the fully relativistic, Oppenheimer-Volkov so- 
lution. This is, because we can choose coordinates such 
that all the terms neglected in the field equations vanish 
identically, and the equations of hydrodynamics (or, in 
this case, hydrostatics) are fully relativistic. Construct- 
ing sequences of equilibrium solutions, we can apply the 
energy turning point method and find the onset of radial 
instability - again without approximation. 

Corotational or irrotational binary stars at large sep- 
arations are very close to being spherically symmetric 
because the two stars interact only through weak tidal 
fields (for example, in Newtonian case, see refs. |p|,^). 
Hence, our formalism can describe the individual stars to 
high accuracy, whether or not the stars are very compact 
and have strong gravitational fields. We therefore expect 
that our approximations are quite adequate. 

In a binary in which the orbital separation a is large 
compared to the stellar radius i?, we can treat the gravi- 
tational effects and tidal deformation due to the compan- 
ion star as a small perturbation. We can then expand the 
gravitational field around the unperturbed, spherical so- 
lution 

a = (0)^ +(2) ae^ +(4) ae^ 

+(5)ae^ +(6) at^ +(7) ae^ + • • •, (2.30) 
V' = (o)V' +(2) V'e^ +(4) V-e** 

+(5)V'e' +(6) V'e' +(7) ^e' + • • (2.31) 
/3' = (i)/3'e+(3) /3^e'+(5) P'e'> 

+ (6)/3^e6 ^^^^ ^^^^ p^^8 + . . (2.32) 

lij ^ % +(2) hije'^ +(4) hije'^ -f (5) /i^j-e^ + • • •, (2.33) 
where we assume 

- - (M/a)i/2 < (i?/a)i/2 _ (2.34) 

and where M is gravitational mass, andjma and (o)V' 
denote the spherical symmetric solutions . Note that 
e denotes the magnitude of gravitational effects from the 
companion star, and hence the expansion in terms of e is 
different from a post-Newtonian expansion (see fl^ .) 

In our formalism, we can calculate (o)Q^ a-nd (o)V' ex- 
actly. Newtonian and first post-Newtonian terms appear- 
ing in („)«, (n)V' a-iid (n)/^* for n =^ arc also taken into ac- 
count consistently for all order in e, although second post- 
Newtonian terms in these and (n)hij are not included. 
Our approximation is therefore appropriate for investi- 
gating the stability of fully general relativistic spherical 
stars due to Newtonian and first post-Newtonian tidal 
effects, which is our goal in this paper. 

Details of our numerical methods, for solving both 
the hydrodynamical equations and the Poisson equations, 



can be found in . We assume symmetry with respect 
to equatorial plane, and solve the equations on a uni- 
form grid of size {2N + 1, 2N + 1, iV -I- 1), covering the 
physical space —L<x,y<L and < z < L where L 
is location of outer boundaries. We use = 50 for the 
spherical symmetric stars, and N = 50, 60, and 75 for 
binary configurations. 

As a numerical check, we monitor the conservation of 
proper mass 

Mp = J p^dV, (2.35) 

total gravitational mass 

M = -2 J S^dV, (2.36) 

and total angular momentum 

J = J {-yJ^ + xJy)ij^dV. (2.37) 

Our difference scheme guarantees conservation of Mp ex- 
actly. Accurate conservation of M and J depends on N, 
and for N ~ 75, the error in one orbital period is ^ 0.05% 
for M and < 0.5% for J, respectively. 



III. DYNAMICAL STABILITY OF SPHERICAL 
STATIC STARS 

In this section we calibrate our code by locating the 
analytically known onset of instability of spherical equi- 
librium stars. 

For initial conditions, we construct sequences of 
equilibria satisfying the Oppenheimer-Volkov equations. 
Note that for these configurations our formalism is ex- 
act and recovers the fully relativistic solutions. In Fig. 1 
we show the proper mass Mp and the (isotropic) coor- 
dinate radius i? as a function of the central density pc 
for a polytrope with F = 1.4 (n = 2.5). We have taken 
advantage of the scale freedom in the problem and cho- 
sen K = 1. Together with c = 1 = G, this assign- 
ment uniquely determines our non-dimensional units for 
length, mass, etc.. For any other value of K, all the re- 
sults can be rescaled trivially (see, e.g., 0). Note that 
for a critical central density pc = Pcrit(— 1.2 x lO^''), the 
mass Mp goes through a maximum Afmax ~ 1.248. This 
maximum marks the onset of radial instability, and sep- 
arates the stable branch {pc < Pcrit) from the unstable 
branch {pc > Pcrit) of the sequence. 

For F = 1.4, the compaction of the maximum mass 
configuration, A/p/i?|crit ~ 0.03, is much less than unity 
(recall Mp/R\„it ^ as F -> 4/3). Nevertheless, the 
mass versus central density equilibrium curve still ex- 
hibits an extremum, and therefore has all the qualitative 
features of any value of F > 4/3 regarding the issue of 
radial stability. Choosing F = 1.4 therefore allows us to 
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FIG. 1. Proper mass Mp (top panel) and isotropic radius 
R (botton panel) as a function of central density pc for rela- 
tivistic, spherical polytopes with F = 1.4. 



study these qualitative features in a regime in which the 
post-Newtonian approximation is very reliable. 

Consider results for five different initial configurations, 
which we denote by A, B, C, D and E (see Fig. 1). We 
construct these initial data with a one-dimensional in- 
tegration of the Oppenheimer-Volkov equations. These 
data are in equilibrium according to the one-dimensional 
finite difference equations. Interpolating these data onto 
the three-dimensional grid of our evolution scheme intro- 
duces a slight perturbation of the equilibrium solution, 
since the truncation error of the three-dimensional finite 
difference equations is different from the one-dimensional 
one. Typically, we find that the pressure of the interpo- 
lated models is slightly larger than the required equi- 
librium value on the three-dimensional grid. We com- 
pensate for that by artificially decreasing the polytropic 
constant K hy a small amount. Note that our finite dif- 
ferencing is convergent, so that for finer grids we need to 
change K by smaller amounts. For the results presented 
here, the initial radius of the star is covered by 30 grid 
points. 

Note that because of the scale freedom in the problem, 
reducing K is equivalent to increasing the mass. This 
is, because iiT"/^, where n is the polytropic index, has 
units of mass (or length) in geometrized units. Consider 
now a model in which we have reduced K by, say, a small 
fraction S <^ 1. We can then rescale this model in order to 
investigate a different K, for example, the original value 
K = 1. Then, the mass of the re-scaled model has to 
increase by a fraction n6/2 to remain in equilibrium. For 
example, reducing K by 0.5 % is, for n = 2.5, equivalent 



FIG. 2. Time evolution of p« at center of stars for models 
A - E. 



to increasing the model's mass by 0.625 %. 

We can now test the stability of our models by vary- 
ing K by small amounts. Small variations in K serve to 
trigger small initial perturbations away from equilibrium 
(e.g. "pressure depletion" ) . Stable models will not change 
their qualitative behavior, whereas unstable models will. 
In Fig. 2 we show the time evolution of the central den- 
sity for the five models. Solid curves are the results 
for K = 0.995, and dotted curves are for K ^ 1. In this 

— 1/2 

section, we plot time in units of Pq , where po is the 
central value of p of the corresponding spherical equilib- 
rium star. 

Obviously, models A and B oscillate stably, for both 
values of K. The period of these oscillations can be com- 
pared with the approximate analytic value [E5l 



tn 



2tt 



^(^-^r3F-4-6.75M 

(5r-6)i?/v R 



-1/2 



(3.1) 



where M is the Newtonian (rest) mass, R is the Newto- 
nian radius of star, and 



pr^dV 



(3.2) 



is the spherical mass moment. Inserting the values pc 
5.6M/R^ and / ~ 0.17MR^ for F = 1.4, we find 



W ^ 5.5p-i/2( 0.2 -6.75^ 



-1/2 



(3.3) 



-1/2 



For model A, M/R ~ 0.017 and hence iosc — ISpc 
which is very close to the value that can be read off in 
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Fig. 2. Obviously, models with a larger compaction M /R 

have, in units of pc ^^"^ , a larger oscillation period tosc, 
which can also be seen in Fig. 2. At the onset of in- 
stability, tosc — oo. From equation (3.3) we therefore 
find that the maximum mass configuration must have a 
compaction M/R ~ 0.029, which is very close to that of 
model C. 

Leaving if = 1 for model C, the star oscillates stably, 
but reducing K by only 0.5 % to 0.995, the central density 
increases monotonically, and the star undergoes gravita- 
tional collapse. This indicates that model C is marginally 
stable against small perturbations, 5K/K < 0.5%, and 
very close to the onset of instability. This is obviously 
true, since its proper mass is only 0.05% smaller than the 
maximum allowed mass Mmax- 

For models D and E, the star monotonically expands or 
collapses, and never oscillates: starting with K = 0.995, 
the star collapses, and iox K — 1, the star expands 
by a large factor Obviously, for if = 1 the star 

should be in equilibrium, and neither expand nor col- 
lapse. The equilibrium is unstable, however, and even 
the smallest truncation error induces a growing pertur- 
bation, which must ultimately lead to gravitational col- 
lapse. Initially, this perturbation may be either an expan- 
sion or a contraction. Since the configuration is gravita- 
tionally bound, the expansion soon has to turn around 
and lead to gravitational collapse [Q. Obviously, we 
would expect that for if > 1 we find an initial expan- 
sion, and for if < 1 an immediate contraction. However, 
due to truncation error, the cutoff between expansion and 
contraction is not precisely at if = 1, but instead at a 
value slightly smaller than unity (if ~ 0.997 when we 
use 30 grid points to cover the star). Again, this value 
approaches unity with increasing grid resolution. This 
behavior establishes that models D and E are unstable 
to radial perturbations, which is, of course, what we ex- 
pected. 

We conclude that we can locate the onset of radial 
instability, and in particular that we can determine the 
maximum allowed mass of neutron stars to very high ac- 
curacy. This will be very important for determining the 
stability properties of binary neutron stars. 



IV. DYNAMICAL STABILITY OF STARS IN 
BINARY SYSTEMS 

In this section we present numerical results on the dy- 
namical stability of stars in binary systems. We always 
assume the two stars to have equal mass, and set up ini- 
tial data so that they are in a circular binary orbit. Note 
that for r = 1.4 the equation of state is sufficiently soft, 
so that there is no innermost stable circular orbit . 
Therefore, we can choose arbitrarily small binary sepa- 
rations, and the orbit of the binary will still be stable. 
We choose very close binaries, in which the separation 
of the surfaces of the two stars is much smaller than the 
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FIG. 3. Proper mass as a function of the maximum den- 
sity for each star in a corotational equilibrium binary (filled 
circles) and for an isolated spherical star (solid line). 



orbital separation {za ~ in the terminology of ref. [Q). 
For these binaries the tidal effects are strongest, and they 
are therefore the most suitable configuration to study the 
stability against gravitational collapse of the individual 
stars. 

We evolve three different classes of binary initial data. 
The first class are corotational binaries. For these con- 
figurations, self-consistent equilibrium initial data can be 
constructed in post-Newtonian approximation i pTf and 
even in general relativity (where the stars are only in 
quasi-equilibrium, see p6|,|7[|). We denote the class of 
post-Newtonian strict equilibrium solutions with a sub- 
script "a". 

In addition to corotational binaries, we would also 
like to study irrotational binaries because they are more 
realistic models for binary neutron stars [^. Self- 
consistent irrotational equilibrium binaries in Newtonian 
gravity have recently been constructed by Uryu and 
Eriguchi |2^, and a relativistic generalization has been 
suggested by Bonazzola, Gourgoulhon and Marck [ pO[ . 
No numerical models are currently available for such 
data in a post-Newtonian approximation. In the absence 
of self-consistent, tidally deformed equilibrium data, we 
therefore take two spherically symmetric stars, put them 
close together, and artificially assign an irrotational ve- 
locity profile which maintains their shape and circular 
orbit approximately. In order to calibrate these irrota- 
tional initial data, which are not strictly in equilibrium, 
we first perform simulations with a second class of corota- 
tional initial data, using spherically symmetric stars, and 
compare these with the self-consistent post-Newtonian 
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equilibrium models which exist in the corotational case. 
For this second class, we assign a uniform velocity 



{-ny,nx,o) 



(4.1) 



to each fluid particle, where is the orbital angular ve- 
locity. We denote these corotational, near-equilibrium 
data with a subscript "b" . 

Finally, the third class of initial data are irrotational, 
near-equilibrium models. Again, we put two spherically 
symmetric stars at a small separation, but now we assign 
an initial velocity 



(4.2) 



{ux,Uy,Uz) = 0,±f2a;o,0 



to each fluid particle |Q . The centers of mass of the two 
stars are located at (±a;o,0,0), where xq > 0. The plus 
sign in (4.2) corresponds to the star at (-|-a;o, 0, 0), and 
vice versa. We denote this third class of initial data with 
a subscript "c" . 

In the above velocity profiles, we determine the angular 
velocity O from Kepler's law 



( Mt 



1/2 



(4.3) 



where Mt — 2Mp is the sum of the proper mass of the 
two spherical stars and a is the coordinate separation 
between their centers of mass. 

In the above velocity profiles, it would be more physical 
to fix instead of Ui. That, however, would involve 
one more iteration in the preparation of the initial data, 
and, in our small compaction cases, would make only a 
negligible difference. We summarize the initial conditions 
for six different models, two in each class, in Table I. 

As in our spherical models, we must vary K slightly 
in order to investigate the stability of the binary mod- 
els. Equilibrium stars in tidal fields are tidally deformed 
and have a slightly smaller central density than spherical 
stars. Using spherical models as initial data for binaries 
therefore overestimates the central density, which causes 
the stars to expand initially. As we will see, this can be 
compensated for by reducing K. 



A. Corotational equilibrium models 

Following Shibata (27) and Baumgarte et al. , we 
construct self-consistent equilibrium initial data, describ- 
ing two corotational binary stars in contact {za — 0). In 
Fig. 3, we plot the proper mass of each star as a func- 
tion of the central density (filled circles), and compare 
these values with those for spherical stars in isolation 
(solid line) . The maximum allowed mass of binary stars 
is slightly larger than that of spherical stars in isolation 
(see the discussion in 0). 

We show results for grid sizes N — 40, 50, 60, and 
75, and find that our code is second order accurate. The 
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FIG. 4. Time evolution of p*max in the corotational equi- 
librium binary models Ba and Da. Solid, dotted and dashed 
lines denote results obtained with A'^ — 75, 60, and 50 grid- 
points, respectively. For model Ba, we set K — 1 initially. 
For model Da, initial values of K are 1, 0.995 and 0.994 for 
iV = 75; 1, 0.992 and 0.99 for N = 60; and 1 and 0.99 for 
N = 50. Curves with numerical labels show values oi K 1. 



proper mass Mp evaluated on a grid with a grid spacing 
d, AIp{d), therefore scales according to 



Mp{d) = Mo + M2(f + Oicf). 



(4.4) 



Taking values for different N (and hence d), we can elim- 
inate the second order error term by Richardson extrap- 
olating to d = 0, which yields the value Mq. This is the 
value that we plotted in Fig. 3. In Table II we summa- 
rize these results by tabulating the masses Mp{d) for the 
different grid resolutions, together with the Richardson 
extrapolated value Mq. 

Comparing Mq with Mp{d) for = 60 and 75, we find 
deviations of ~ 0.6% and ^ 0.4%, respectively. This is a 
lower limit on the truncation error that we have to expect 
in the subsequent evolution. 

Note that in Fig. 3, we plot the mass versus maximum 
density for a sequence of constant separation (namely for 
contact binaries). In this graph, the onset of instabil- 
ity need not coincide with the maximum mass configura- 
tion. Instead, the onset of instability can be located by 
constructing sequences of constant angular momentum 
(see Baumgarte et al. jT^). However, we expect that the 
onset of instability is very close to the maximum mass 
configuration. We therefore present results for two mod- 
els, one with a maximum density slightly less than the 
maximum mass model (denoted Ba, Pmax = 5 x 10~^), 
and one with a maximum density slightly larger (Da, 
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FIG. 5. Snapshots of density contour lines and the velocity 
flow {v^jV"") in the equatorial plane for model Ba. Contour 
lines are drawn for pt/p, o = lO"" "^-' , where p* o denotes 
the maximum value at t = 0, for j = 0, 1, 2 ■ ■ ■ 10. Vectors 
indicate the local velocity field. Time is shown in units of 
orbital period P. See Table I for the relation between P and 

Pmax,0 . 

Pmax = 1-5 X lO^"'). For both models, the orbital period 

— 1/2 

is 50 in units of p^^J^ q where Pmax,o is maximum value 
oi p at t = 0. 

In Fig. 4, we show the time evolution of the maximum 
value of (/5*max) for models Ba and Da. We show 
results for three different grid resolutions, N — 50, 60, 
and 75, where we have kept the location L of the outer 
boundary constant. We picked L such that each star is 
covered by ~ N/2 — 5 grid points (see Table I). 

For model Da, we also picked several different values of 
K, and marked all simulations with K ^ 1 accordingly 
in Fig. 4. Depending on K, these models either collapse 
or expand, but never oscillate stably. This indicates that 
they are dynamically unstable, as we have expected. As 
in the discussion of unstable spherical models, we would 
expect the cutoff between initial expansion and contrac- 
tion to be at if = 1 if we had arbitrary accuracy. This is 
not the case, but we again find that increasing the grid 
resolution makes this cutoff approach unity (cutoff value 
of K is less than 0.99 for N = 50, ~ 0.992 for N = 60, 
and ~ 0.995 for N = 75). 

For model Ba, we only show results for K = 1. Obvi- 
ously, this configuration oscillates stably. The oscillations 
are due to a slight inconsistency between the initial data 
and the evolution scheme. The amplitude of these oscilla- 
tions decreases with increasing grid number, which shows 
that our method is convergent. Note that Mq of this 



FIG. 6. Same as Fig. 5, except for model Da at t = 0. For 
this sequence we set K = 0.994 initially. 

configuration is 1.250, which is marginally larger than 
the maximum allowed rest mass of an isolated, spherical 
star, Afinax = 1.248. We therefore conclude that all stars 
that are stable in isolation are also stable in a corota- 
tional binary. We expect that the reverse is not true: a 
star in a close binary can support more mass than an iso- 
lated, spherical star. Note that our results are different 
from those of WMM, who found that stars with masses 
of as much as 10 % or more below the maximum allowed 
rest mass in isolation were destabilized in a close binary. 
At this level, such an effect would be discerned easily by 
our code, but it is not present. 

In Figs. 5 and 6 , we show contour lines of and the 
velocity field of {v^,v'^) in the equatorial plane for mod- 
els Ba and Da, using a grid resolution of iV = 75. For 
model Da we set K — 0.994. Note that for our adopted 
soft equation of state, the stars are very centrally con- 
densed. The tidal field mostly deforms the very low den- 
sity envelope, which is hardly visible in these plots. The 
two envelopes are pulled towards the binary companion, 
and in our contact cases, touch at the origin. The high 
density cores, however, are hardly deformed by the tidal 
field, and are still fairly far separated. 

For model Ba, we show contours at t = 0.0, 27.0, and 

— 1/2 

55.5, all in units of p^^/^ g, which corresponds to the ini- 
tial condition, a little more than half an orbital period, 
and a little over an orbital period. It is obvious from the 
graphs that the two stars stably orbit each other. Note 
also that the velocity field remains close to being corota- 
tional, and that we do not see any evidence of vorticity 
features which have been reported in ref. |8| . 

For model Da we show contours at t = 0.0 9.55 and 
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FIG. 7. Time evolution of p,max in the corotational, 
near-equilibrium binary models Bb and Db. Initially, we set 
K = 0.98 (dotted line), K = 0.975 (solid line) and K = 0.97 
(dashed line) for model Bb, and K = 0.98 (dotted line) and 
K = 0.975 (solid line) for model Db. For models Bb and Db, 
the orbital period is 



39 and mp-lil. 



— 1/2 

19.1p,jjjj^ Q, which is a httle less than half an orbital pe- 
riod. It can be seen very clearly how the star contracts 
and collapses. 



B. Corotational near-equilibrium models 

We now present numerical results for our corotational 
near-equilibrium models. We do this to calibrate our 
code, and to show that our near-equilibrium models are 
goo d ap proximations to self-consistent equilibrium mod- 
els p2|. This justifies studying such near-equilibrium 
models for irrotational binaries. 

For all the simulations discussed in this and the follow- 
ing subsection, we use a numerical grid with = 60 grid 
points. We adjust the outer boundary L such that the 
radius of each star is covered with ~ 25 grid points. For 
the models Bb and Be we used L — 125 and for Db and 
Dc L = 96 (see Table I). We construct spherical mod- 
els, and placed them on the grid such that their centers 
of mass is located at {x,y,z) = (±L/2, 0,0). Note that 
these stars are not in contact. However, the stars are not 
tidally deformed, and therefore separation of the center 
of masses of the two stars is slightly smaller than for 
the self-consistent equilibrium models. The orbital pe- 



riod for these binaries is 
respectively. 



39 and 46 in units of p 
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FIG. 8. Same as Fig. 5, but for model Bb. For this sequence 
we set K = 0.975 initially. 



In Fig. 7, we show the time evolution of p^max for mod- 
els Bb and Db. Dotted, solid and dashed lines denote 
results for K = 0.98, 0.975 and 0.97, respectively. As 
before, model Db either expands or contracts, but never 
oscillates stably. We again conclude that Db is dynam- 
ically unstable. Model Bb, on the other hand, exhibits 
stable oscillation for several different values of K, and we 
conclude that it is dynamically stable. 

It is interesting to note that even with a reduced value 
of K = 0.98, model Bb initially expands, albeit stably. 
Only reducing the pressure to a value smaller than that 
{K ^ 0.975), the configuration is roughly in equilibrium. 
This can be understood quite easily, because the spheri- 
cal star is not a self-consistent equilibrium solution. The 
tidal field tends to deform the star, which reduces the 
central density. Therefore, our initial data have too high 
a central density for their mass, and the star starts ex- 
panding. This can be compensated for by reducing K and 
hence the pressure. As we have argued before, reducing 
K is equivalent to increasing the mass, and indeed this 
is another way of understanding why a star in a binary 
can support more mass than a star in isolation. 

In Fig. 8, we show contour lines of and the velocity 
field {v^ jV"^) in the equatorial plane for model Bb, where 
we have set K = 0.975. We show the configuration at 

— 1/2 

t = 0.0, 26.5 and 53.0pj^3^^ q. Comparing these plots with 
Fig. 5, one can see that the center of masses of the two 
stars are closer. During the evolution, the stars loose 
their spherical shape and adjust to the tidal field. How- 
ever, all the qualitative features are very similar to the 
self-consistent equilibrium simulations. In particular, the 
velocity field remains nearly corotational, and we do not 
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FIG. 9. Time evolution of p*max in the irrotational, 
near-equilibrium binary models Be and Dc. Initial values of 
K are 0.99, 0.985 and 0.98 for model Be, and 0.98 and 0.985 
for model Dc. Solid, dotted and dashed lines denote results 
for K = 0.98, 0.985 and 0.99, respectively. The orbital period 

1/2 

of models Be and De is ~ 39, and 46^^^^^^^ q, again. 



see any evidence of the double vorticity fields as reported 
by WMM. 

This test suggest that the spherical near-equilibrium 
models are very good approximations to self-consistent 
equilibrium configurations for the models and issues we 
are investigating here (but see footnote P3l). 



C. Irrotational near-equilibrium models 

Viscosities in neutron stars are expected to be too 
small to bring neutron star binaries into corotation on the 
timescale of their evolution [ 28[ . Numerical models exist 
only in Newtonian gravity | 29f , |30| or in the Newtonian 
and post-Newtonian ellipsoidal approximation p2| , pl| . 
We therefore adopt the spherical near-equilibrium ap- 
proximation to construct initial models, which we have 
calibrated and found to be adequate in section IV. B for 
corotational cases. 

We again vary K and choose K = 0.99, 0.985 and 0.98 
for model Be, and K = 0.985 and 0.98 for model Dc. In 
Fig. 9, we show time evolution of /0*max for models Be 
and Dc. Dashed, dotted and solid lines denote results for 
K = 0.99, 0.985 and 0.98, respectively. As in the coro- 
tational cases, model Dc cannot be held stably, whereas 
model Be oscillates for all these choices of K. Therefore, 
we again conclude that model Be is dynamically stable, 
whereas Dc is not. 



FIG. 10. Same as Fig. 5, but for model Be. For this se- 
quence, we set K = 0.99 initially. 



It is interesting to note that for these irrotational mod- 
els, we had to reduce if by a somewhat smaller amount 
to minimize the amplitude of oscillations in model Be 
than for the corotational model Bb {K ^ 0.985 here, 
and K ~ 0.975 for model Bb). This can be understood 
very easily, because in corotational models, the individ- 
ual stars are spinning, and are therefore stabilized and 
deformed by both the tidal field and the their own spin. 
In irrotational binaries, the stars have almost no spin 
(with respect to distant inertial observers), and are de- 
formed only by the tidal field. Therefore, putting a star 
into an irrotational binary will reduce the central den- 
sity by less than putting the same star into a corota- 
tional binary. As a consequence, we have to reduce K by 
a smaller amount to compensate. Applying our scaling 
argument, this result means that a irrotational binary 
can support less mass than a corotational binary, but 
still more than a spherical star in isolation. This result 
is corrobated by the post-Newtonian ellipsoidal models 
constructed in ||ll[| . 

In Fig. 10, we show contour lines of and the veloc- 
ity field (d^, v^) in the equatorial plane for model Be and 
K — 0.99 initially. We show contours at t — 0.0, 26.5, 

— 1/2 

and 53.0pjjjg^^ Q. For the bulk of the matter at the core 
of the stars, the velocity field remains approximately ir- 
rotational, and the stars stably orbit each other. We 
conclude that there is no qualitative difference between 
corotational and irrotational binaries as far as their radial 
stability properties are concerned. 
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V. SUMMARY 
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We perform post-Newtonian, dynamical simulations of 
close binaries in circular orbit. In particular, we study 
the stability of the individual stars against gravitational 
collapse in both corotational and irrotational systems 
containing stars of equal mass. 

We have chosen a soft, adiabatic equation of state with 
r = 1.4, for which there is no innermost stable circular 
orbit, so that the binary orbit is stable even when the 
stars are in contact, and for which the onset of instability 
for a spherical star in isolation occurs at a very small 
value of the compaction M/R. We can therefore study 
the individual stars' stability properties in near contact 
binaries, for which the tidal effects are strongest, and in a 
regime in which a post-Newtonian approximation is very 
accurate. 

We do not find any crushing effect as reported by 
WMM [^. In contrast, the maximum density in both 
corotational and irrotational binaries is smaller than that 
of spherical stars in isolation. We find that stars in bina- 
ries can support more mass than in isolation. Moreover, 
all stars that are stable against radial perturbations in 
isolation, will also be dynamically stable when put into 
a binary. 

All these results are in complete agreement with, for 
example, the findings of Baumgarte et al. ^6 17\, Flana- 
gan iQ, and Thorne ||l^. For the most part, their dis- 
cussions rigorously address secular stability only. Several 
different arguments can be invoked to suggest that secu- 
larly stable binaries are also dynamically stable, but this 
is strictly proven only in Newtonian theory (see also ) . 
Our dynamical calculations reported in this paper are 
the first to directly confirm dynamical stability, at least 
within our post-Newtonian approximation. 

We compare, in a near-equilibrium approximation, 
corotational and irrotational binary models. As ex- 
pected, stars in corotational binaries can support slightly 
more mass than in irrotational binaries, but apart from 
these small differences we do not find any qualitative 
difference in their radial stability properties. A more 
rigorous treatment will require the construction of post- 
Newtonian, irrotational equilibrium binary models for 
initial data. 

Since our computations have been performed in nondi- 
mensional units, our results apply not only to neutron 
star binaries, but also to binaries of white dwarfs and su- 
permassive stars. In fact, the equations of state of mas- 
sive white dwarfs (ideal degenerate, extremely relativistic 
electrons) and supermassive stars (radiation 3> thermal 
pressure) are closely approximated by the value F = 1.4 
that we have adopted. These binaries may be impor- 
tant low frequency gravitational wave sources for future 
space-based gravitational waved detectors, like LISA. 
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Table I. Initial conditions of binary models. We tabulate the maximum density Pmax.o/10 ^, the radius i?oo of 

— 1/ 2 

a spherical star of the same rest mass in isolation, the orbital period P in units of p^^^ q, the nature of the initial 
velocity field and matter profile, and the location of the outer boundary L (see text). All quantities are shown in 
units oi c = G = K = 1. 



Model 


Pmax,o/10 ^ 




f X -^Pmax.O 


Velocity Field 


Matter Profile 


L 


Ba 


5 


53 


48.6 


corotation 


equilibrium 


141 


Da 


15 


38 


50.2 


corotation 


equilibrium 


97.5 


Bb 


5 


53 


39.4 


corotation 


spherical 


125 


Db 


15 


38 


45.9 


corotation 


spherical 


96 


Be 


5 


53 


39.4 


irrotation 


spherical 


125 


Dc 


15 


38 


45.9 


irrotation 


spherical 


96 



Table II. Proper mass Mp of each star in a corotational equilibrium binary for different grid resolutions. We tabulate 
Mp versus central density Pc for N = 40, 50, 60, and 75. We also list the Richardson extrapolated value Mq as well 
as the proper mass of the spherical model in isolation with the same central density. 





2 


5 


10 


15 


N = m 


1.1954 


1.2347 


1.2469 


1.2451 


iV = 50 


1.2004 


1.2402 


1.2526 


1.2509 


iV = 60 


1.2031 


1.2431 


1.2556 


1.2539 


AT = 75 


1.2051 


1.2455 


1.2581 


1.2565 


Mo 


1.209 


1.250 


1.263 


1.261 


Spherical 


1.1933 


1.2344 


1.2482 


1.2476 
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